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Abstract 

Several cellular automata (CA) models have been developed to sim- 
ulate self-organization of multiple levels of structures. However, they do 
not obey microscopic reversibility and conservation laws. In this pa- 
per, we describe the construction of a reversible lattice molecular au- 
tomata (RLMA) model, which simulates molecular interaction and self- 
organization of higher-order structures. The model's strict reversibility 
entails physically relevant conservation laws, and thus opens a way to pre- 
cise application and validation of the methods from statistical physics in 
studying the necessary conditions for such multiple levels of self-organization. 

Keywords: Reversible Cellular Automata; Molecular Aggregation; Self- 
organization; Artificial Chemistry 

1 Introduction 

Possessing and utilizing multiple levels of self-organized structures — sometimes 
addressed as dynamical hierarchies[l, 2] — is a characteristic feature of biologi- 
cal systems. Cellular automata (CA) and similar discrete paradigms have been 
effective in modeling such dynamical self-organization hierarchies. In the con- 
text of molecular aggregation, lattice molecular automata (LMA) simulates self- 
organization of water (polar solvent), monomers, and polymers into clusters and 
higher-order structures such as miccllcs[3, 4, 5], and similar models have been 
developed to simulate organization of compartment structure and proto-cell-like 
self-reproduction [6, 7]. 

However, these models do not obey microscopic reversibility and conser- 
vation laws, and therefore, the possibility and stability of the self-organized 
structures in these models are, to some extent, implied in their irreversible time 
evolution rules. Under the laws of physics, stable persistence of an organized 
structure requires effective utilization of limited resources and smooth disposal 
of generated entropy. Therefore, the constraint of reversibility should not be 
omitted in studying the necessary conditions for stable structures, using, for 
example, the canonical methods of statistical physics. 
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In this paper, we describe the construction of reversible lattice molecular 
automata (RLMA), which simulates self-organization of water, monomers, and 
polymers with a strictly reversible dynamics and physically appropriate conser- 
vation laws. Although several reversible CA models have been proposed to sim- 
ulate self-organization processes[8, 9], our model can simulate self-organization 
of structures with mobility, which will be essential for realizing higher-order 
structures and higher functionality such as autonomy. 

The rest of the paper is organized as follows. A formal definition of CA is 
provided and useful techniques in constructing reversible CA are briefly reviewed 
in Section 2. Features of the original LMA model and its relation to other 
models are also briefly reviewed there. The construction of our RLMA model 
is described in Section 3, along with the conservation laws derived from the 
reversible dynamics. Some simulations of monomers and polymers in polar 
solvent are presented in Section 4. Finally, our conclusion is drawn in Section 
5. Appendix A presents an alternative approach for implementing reversible 
molecular rotation. 

2 Preliminaries 

2.1 Formalization of CA 

On a d-dimcnsional spatial lattice Z d , each site (cell) i G Z d is assigned with 
its local state oji G A. The finite set A of local states is called an alphabet. A 
specification of local states over the whole space oj = (uji) ie %d G Q = A % is 
called a global state or configuration. 

The dynamics of a CA is given by the local transition map ip as 

where the neighbor function AT : i i— ► . . . ,j N ) defines interaction range for 
each site i. By applying the local map ip over the lattice, the global transition 
map 

Lu t+1 = <P(w*) (2) 

from a configuration at t to the one at t + 1 is derived. Although the applica- 
tion of ip over the space is synchronous in simple CA, making it asynchronous 
can be effective in satisfying reversibility and other constraints, as shown later. 
Furthermore, in more complicated CA, the local map (p consists not of a single 
map but of several maps (sub-steps), and the local states u>i G A also have inner 
structures such as "partitions" or "layers." 

When the global transition map ^ is bijective, that is, when for any config- 
uration uj t+l its pre-image u>* is unique, the CA is reversible (or invertible). 

2.2 Construction of reversible CA 

Reversibility entails conservation of information — differences in states cannot 
just appear from or vanish into nowhere. Hence the manner in which to prevent 
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information loss is crucial in constructing reversible CA. 

Since the many-to-one local transition map if in Eq. (1) obviously loses infor- 
mation by itself (Fig. 1(a)), the loss should be prevented by well-counterbalanced 
distribution among the interacting cells. However, designing such maps is far 
from trivial. Indeed, judging the reversibility of a global transition map <P, given 
its corresponding local map if, is a difficult task in itself. 

An easier method to construct reversible CA is by adopting a permutation 
(reversible transformation on a finite set) ip, 

iP:A B ^ A B , (3) 

instead of a many-to-one mapping if, as a constituent of the transition rule 
(Fig. 1(b)). Here, B denotes a "block" of cells under the permutation. In 
Partitioning CA (or block CA) [10], for example, both the reversibility and global 
transmission of information are satisfied by combining the permutation and 
alternation of different partitioning schemes of a given space into cell-blocks. 
The permutation (3) can be generalized into a conditional permutation as 

ij c :A s ^A s , ceA c . (4) 

Out of the set C + S of cells that are subject to the mapping, the states of cells 
in C work as "control signals," which determine a permutation for the states of 
S, and reappear unchanged as outputs (Fig. 1(c)). The Fredkin gate and the 
Toffoli gate arc well-known examples of conditional permutations. 



(a) (b) (c) 




Figure 1: Mappings to constitute transition: (a) many-to-one mapping, (b) 
permutation, and (c) conditional permutation. 

To prevent information loss, the outputs of the (conditional) permutations 
should be reused as inputs or conditional signals, and to achieve this, one needs 
to implement certain techniques such as dividing a time step into several sub- 
steps, and arranging the permutations sparsely and asynchronously in space- 
time[ll]. 

For synchronous information transmission, one can also use global shifts, 
which uniformly displace some partitions or layers of local states. The transla- 
tional movement of free particles can be effectively modeled by the shifts. In 
lattice gas automata (LGA)[12, 13], for example, shifts are utilized to express 
the translation of gas particles, in combination with permutations that repre- 
sent the collisions of the particles. Partitioned CA [14] is virtually equivalent to 
the LGA. 
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2.3 LMA model and other models 



Various models have been proposed and used to simulate molecular self-organization 
processes. On the one hand, the molecular dynamics (MD) method models 
molecules as particles with appropriate interaction potentials, and solves their 
equations of motion in continuous space[15]. While the MD enables micro- 
scopically detailed description of the dynamics, size of the simulated system is 
restricted by the available computer resources. On the other hand, lattice-type 
models have been successful on simulating macroscopic behavior of phase sepa- 
ration and aggregation processes. Especially the Larson model[16, 17, 18] and 
its variants[19, 20, 21, 22] are widely used and many results are reported. In the 
traditional Larson model simulating ternary mixture of water, oil (hydropho- 
bic monomers), and surfactant (amphiphihe polymers), water and hydrophobic 
monomers are represented by a set of up and down (+1,-1) spins, respectively, 
and polymers are represented by strings of spins. Monte Carlo method is used 
for update and the ferromagnetic interaction between the spins induces phase 
separation, micelle formation, etc. 

The original LMA model bridges the gap between the MD method and the 
Larson- type models [5]: While realized in a discrete manner and thus keeping the 
efficiency of the lattice setting, it includes some microscopic molecular details, 
such as hydrodynamics conserving momenta in the molecular collision, direc- 
tions of polar molecules and accompanied anisotropy of molecular potential 
energy. A distinguishing feature of the LMA model is the equienergetic inter- 
action for the pairs water-hydrophobic monomer and hydrophobic monomer- 
hydrophobic monomer, following experimental data on enthalpy exchanges in 
mixtures [23]. This setting is in contrast to the Larson models, which define 
positive enthalpic gains for oil-oil interaction but not for water-oil interaction. 
Consequently, in the LMA model phase separation is realized via entropy- driven 
hydrophobic effect, and not enthalpy-driven as in the Larson-type models. 

Although update rule of the LMA model partially keeps the conservation 
laws, its dynamics is not microscopically reversible (refer to section IV. B and 
V of Ref. [3] for example to see the total energy is conserved in the mean but 
not strictly and explicitly). Therefore, utilizing the techniques introduced in 
section 2.2, we construct our RLMA model in the next section. 

3 RLMA Model 

3.1 Space 

We formalize the RLMA model on the two-dimensional triangular lattice (Fig. 2(a), 
(b)) as in the literature[3, 4], although generalization to other lattice struc- 
tures and to higher dimensions will be straightforward. We use the variable 
I G {+1, +2, +3, —1, —2, —3} = L to denote the principal directions, and to 
denote cell i's nearest neighbor in direction I, as shown in Fig. 2(c). L corre- 
sponds to {0, 7r/6, . . . , 57r/6} in the equilateral triangular lattice with a proper 
coordinate system (Fig. 2(c)), and on L, wc define a cyclic permutation A of 
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length 6, 



A = 



-1 +2 +3 -1 -2 -3 
-2 +3 -1 -2 -3 +1 



(5) 



which corresponds to +7r/6 rotation operator for the principal directions. 

(c) 



(a) , (b) 
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Figure 2: Two-dimensional triangular lattice: (a) structure, (b) corresponding 
cells, and (c) principal directions and nearest neighbors. 



3.2 Local states 

Each local state has the layers (internal data structure) shown in Table 1. 



Table 1: Structure of local states. 



Layer name and variable 


Values 


Molecular type mti 


Water (W), hydrophilic monomer (I), 
hydrophobic monomer (0), or vacuum 
(V) 


Molecular orientation moj 


moi G L for polar molecules, moi = 
null otherwise 


Translational kinetic energy (tka t i)ieL 


tkeij <E {0, 1} for molecules, while non- 
zero values in the opposing directions 
are forbidden 


Rotational kinetic energy rkei 


{ — 1, 0, +1} (polar: {±1}, non-polar: 0) 


Molecular bonds (mbi t i)i^L 


Up to two bonds for hydrophilic or hy- 
drophobic monomers 


Heat particles {hi^i^L 


hi t i € {0, . . . , H max } for each I E L 


Preferential direction pdi 


pdi e L 



For each cell i £ Z 2 , molecular type mti takes one of three types of molecules — 
water (W), hydrophilic monomer (I), hydrophobic monomer (O), or vacuum 
(empty; V). For example, one can consider the hydrophilic monomer to be 
acetic acid and the hydrophobic monomer to be methane. A site can contain at 
most one molecule; this constraint corresponds to excluded volume. 
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Figure 3: Molecular types. Upper row: schematic illustrations. Lower row: 
representations in visualization of the simulation results; bars indicate the ori- 
entations of polar molecules. 



Water and hydrophilic monomers are polar molecules; therefore, they have 
molecule orientation moi € L. (For hydrophobic monomers and vacuum, 
moi = null.) We define that, for water in orientation mo, the same direction 
represents negative polarization (corresponding to one oxygen) and A ±2 (mo) 
represents positive polarization (corresponding to two hydrogens), and for a 
hydrophilic monomer in orientation mo, A ±1 (mo) represents negative polariza- 
tion (corresponding to O or OH) (See Fig. 3). Molecular orientation affects the 
strength of potential energy induced by several kinds of molecular interaction 
(see section 3.3). 

The sites occupied by molecules have translational kinetic energy (TKE) 
tkei.i G {0, 1} in every principal direction I G L, although non-zero energy 
values in opposite directions on the same line are forbidden {tke^i + tke^-i < 1). 
Hence, there are 27 possible TKE states for a molecule (Fig. 4). 




Figure 4: Possible states of translational kinetic energy for a molecule. 



Molecules can have rotational kinetic energy (RKE) rkei, which allows the 
rotation of the polar molecules to be reversible (see section 3.4.4). For proper 
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update by the rotation rule given here, we confine the value of rka to {+1, —1} 
for polar molecules, and to zero for non-polar molecules or vacuum. (For an 
alternative setting, see Appendix A.) 

Hydrophilic and Hydrophobic monomers can have molecular bonds with 
neighboring monomers. We define that 

, f 1 if two monomers at i and l(i) are bonded, ,„■. 
mh ' 1 = { otherwise. (6) 

(Thus, mbij = m&j«) Polymers can be composed as a group of monomers 
linked by the bonds, as shown in Fig. 5. In the current study, we suppose that 
for each monomer to have the bonds in at most two directions, X^zei m ^i,i — 2; 
thus, the polymers are one-dimensional. One can consider the polymers to be 
fatty acids. 
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Figure 5: Schematic representation of a polymer. 



For the above layers, which are related to molecules, we also overlay the heat 
particles layer on each cell. The heat particle variable h+j can take values of 
{0, . . . , iifmax} independently for every direction I G L. 

Finally, we append the preferential direction pdi G L for each cell i. In the 
transition rule given below, the preferential direction works as a "fluctuation" to 
break irreversibility-inducing symmetry. The parity of the preferential direction, 
defined by 

-1 if pd i G {+1, +3, —2}, 
-1 if pdi G {+2,-1,-3}, 



parity (pdi ) = 



(7) 



is also utilized in the transition rule. 

Molecular type, orientation, TKE and, molecular bonds (or variables equiv- 
alent to them) are included in the original LMA[4, 5]. On the other hand, 
RKE, heat particles, and preferential direction are introduced in this model to 
implement reversibility in a physically appropriate manner. 
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3.3 Potential energy 

Every molecule interacts with its nearest neighboring molecules; 1 therefore, it 
has potential energies for each of the six principal directions. In calculating 
potential energy, we consider only pairwise interactions, and let V" 1 ' 1 ^ denote 
potential energy arising from the interaction between molecules at i and l{i). 
The molecular interaction is divided into three classes: 2 

• Electrostatic interactions between permanent multipoles, which take place 
when the polarized directions of the two polar molecules face each other. 
Let V'pcrm-pcrm represent the potential energy contribution from this class 
of interactions. 

• Induction-based interactions between a permanent multipole and an in- 
duced multipolc, which take place when a polarized direction of one molecule 
faces an originally non-polarized direction of another. Let Vp erm _i n( j rep- 
resent the potential energy contribution from this class of interactions. 

• London dispersion interactions between induced multipoles, which take 
place when the surfaces of two non-polar molecules face each other. Let 
Vind-ind represent the potential energy contribution from this class of in- 
teractions. 

Then, the total potential energy in the system is calculated as 

v* t , = I W v iMi) = - W (v im +v im +v im \ 

'total g / j / j 2 L < ' j \ perm— perm ' " perm— ind 1 ind— indy ' 

i leL i leL 

(8) 

For the full specification of the potential terms in our model, the integer 
parameters listed in Table 2 must be given. 

3.4 Transition rule 

In the original LMA model, each unit-time update consists of the following 
sub-steps [4]: 

1. propagation of the molecular type and redistribution of kinetic energies, 

2. construction of type-specific force fields, 

3. calculation of potential energies, 

4. calculation of the most proper move direction, 

5. readjustment of bonds in polymers according to the move direction, and 

6. movement of the molecule and clearing of the old lattice position. 

1 Although wider range of interaction can also be modeled, it requires larger number of site 
groups and more complicated update schemes (see section 3.4). 

2 We omit cooperativity effects because of their minor influence on the simulation results. 
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Table 2: Parameters of potential energy. Potentials of other neighboring direc- 
tions of molecular pairs are set to 0. 



Class 


Potential 


Applied cases 


■"perm — perm 


V\VH-WH 

Vwo-wo 

VwH-WO 


Where positively polarized directions (Hs) of two water 
molecules (Ws) face each other 

Where negatively polarized directions (Os) of two Ws 
face each other 

Where an H and an of two Ws face each other 


V\VH-IP 
V\VO-IP 

Vip-ip 


Where an H of a W faces a negatively polarized direc- 
tion (0 or OH) of a hydrophilic monomer (I) 
Where an of a W faces a negatively polarized direc- 
tion of an I 

Where two negatively polarized directions of two Is face 
each other 




VwH-WN 
V\VH-IN 


Where an H and a non-polarized direction of two Ws 
face each other 

Where an H of a W faces a non-polarized direction of 
an I 

Where an H of a W faces any one of directions of a 
hydrophobic monomer (0) 

Where an and a non-polarized direction of two Ws 
face each other 

Where an of a W faces a non-polarized direction of 
an I 

Where an of a W faces any one of directions of an 
Where a negatively polarized direction of an I faces a 
non-polarized direction of a W 

Where a negatively polarized direction and a non- 
polarized direction of two Is face each other 
Where a negatively polarized direction of an I faces any 
one of directions of 


Vperm— hid 


VwH-O 

V\VO-WN 

VwO-IN 




Vwo-o 

^IP-WN 
VlP-IN 

Vip-o 


Vind— ind 


Vo-o 


Where any directions of two Os face each other 



Although stated otherwise in Ref. [5] , many of these sub-steps are irreversible 
in actuality, involving erasure and duplication of information. To realize re- 
versibility, therefore, we reconstruct the sub-steps and create new ones, utilizing 
the techniques introduced in section 2.2. 

3.4.1 Molecular translation, collision, and excluded volume 

In the LMA, for each molecule, the most proper move direction is calculated 
based on its TKEs and potentials, and the molecule moves to the direction if 
the movement satisfies the constraints of excluded volume and molecular bond 
maintenance. This rule causes situations whose prc-images are not unique (e.g., 
a molecule at a site might have come from one of the neighboring sites according 
to the most proper move direction, or might have been at the same site a unit 
time ago because of the constraints), and thus it is irreversible. 

To satisfy the constraints of excluded volume and molecular bond mainte- 
nance, and to realize reversibility at the same time, we introduce site groups. 
Sites in each group should be scattered uniformly and sparsely enough (to pre- 
vent interference of the pairwise interactions defined below, the sites in each 
group should be separated by at least four times the unit distance). We deter- 
mine the group to which site i belongs at time t by the following map 

g(i, t) = {4(i x mod 4) + (i y mod 4) + t} mod 16, (9) 

and let G = {0, 1, . . . , 15} denote the range of g. Here (i x , i y ) are the coordinates 
of site i given by the axis in Fig. 2(c), and each site is assigned to a group, as 
shown in Fig. 6. 

(7) © © © ® © 

© ® @ (D © 

©©©©©© 

\ / \ / \/\ ' 

(0) (4) (&) © (0) 

x / \ / \ / \ / \ / \ 
® ® ® ® © © 

\ / \ / \ / \ / \ / 

— ©> — w — ® @ 




Figure 6: Site groups for interleaved interaction. Sites in the group indexed "6" 
arc shaded. 



Using the site groups and the preferential directions, molecular translation 
and collision are performed in an interleaved manner using the scheme shown 
in Fig.7. 
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begin 

for g in site groups G do 
for k in (0,1,..., 5) do 

for every i in a group g do 

3 ~pdi{i); 

if A k (pdj) — pdi then 

^ mtc (cjj , ) ; / /pairwise interaction between i and j. 

end 



Figure 7: Interleaved paired site interaction scheme. Note Eq. (5) for the defi- 
nition of A. 



For the interaction between paired neighboring sites i and j = pd^i) in 
Fig. 7, we define a composite conditional permutation i/j m t c , which represents 
molecular translation and collision. 

When both sites are vacuum, no interaction takes place: 



rati = mtj = V => Vmtc(^i, Uj) 



{uJi,LUj) (identity). 



(10) 



When i is occupied by a molecule and j is vacuum, the molecule at i moves 
to j if the TKE in this direction pd i is positive, or the molecule's TKE in 
this direction is inverted if the TKE in the opposite direction is positive, or no 
change takes place if the molecule's TKEs are zero in both directions: 



rat, 



7^ V, ratj 
( ( 



= V = 
rati 
raoi 
{tke it i)i 
rkti 
V {mbi,i)i J 



( mtj \ 

moj 

rkej 
V (mbj,i)i J 



ipmtcitkei.pd^^, tkei^pd^ 



/ mtj \ 
moj 

( tke 3,l)l 

rkej 

\ \ ( mb jd)i J 
if (tke 



mti 
mOi 

rkei 
V {™>bi,l)l J 

1*0 = (i.o) 



(tkCi^ — pdj^i tkGi^pd^ 

if {tkei tPd ., tkei-p d .) = (0, 1), 

if {tkei.pd.,tkei._p d .) = (0,0). 

(11) 

(Here, internal layers that arc unaffected by ip m tc are omitted.) Conversely, 
when i is vacuum and j is occupied by a molecule, the molecule's move direction 
pdi is replaced by — pc? i in the above equation. Note that if the preferential 
directions arc uniform over all of the sites, a free molecule without molecular 
interaction always maintains the directions of its TKEs after one full unit-time 
update, because whenever an inversion of TKE takes place, it is canceled by 
another inversion in the update scheme of Fig. 7. This does not hold, however, 
if the preferential directions are not uniform. This issue will be addressed again 
in section 3.6. 
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When both i and j are occupied by molecules, the two molecules exchange 
TKEs as in an elastic collision 



mti^V A mtj ^ V ^ 
3.4.2 Maintenance of molecular bonds 

As mentioned earlier, polymers are composed as chains of monomers linked by 
molecular bonds (Fig. 5). Integrated and coherent motion of such a multi-site 
structure ("solid body") is difficult to model using CA. A possible approach is to 
express the structure's motion states by its deforming shape[24, 25]. Although 
this method can replicate many aspects of Hamiltonian mechanics as well as the 
structure's integrity, it makes it difficult to formalize proper interaction between 
such a structure and the single-site particles (molecules) whose motion states 
are expressed as their internal states. In the RLMA, the maintenance of bonds 
is ensured by using the bond information as another conditional signal for the 
molecular translational permutation (11). 

First, it is checked if the molecular bonds (provided they exist) are main- 
tained when the molecule at i moves to the vacuum site j = pd^i). Bonds 
with molecules at A zizl (pd i )(i) are not destroyed by the movement because the 
bonded molecules remain neighbors (the bond directions after movement be- 
come A ±2 (pd i ); Fig. 8(a)). If the molecule has bonds in other directions, these 
bonds are destroyed by the movement (Fig. 8(b)). Therefore, we append another 
condition to the permutation (11): 

fnti 7^ V, mtj = V 

Apply (11) followed by bond readjustment -^br if «i&i.A± 2 (p(2 ) = m bi,-pdi = 0, 
i>mtc(ui, uij) = {u)i, ujj) otherwise. 

(13) 

The readjustment of bonds V'br takes place not only at the moved molecule's 
new position j but also at A ±1 (pd i )(i) if the bonds exist: 



V'br 



m 



bj,A- 1 {pd i ), m bj,A- 2 {pd i ) 
\ m bA- 1 (pd i )(i),A+ 2 (pd i )7 m bA- 1 (pd i )(i),A + 1 (pd i ) I 



0, m& A + i( pd .) (i ) ,A- 2 (pd 

0, mbj^-i( pd .) 

\ 0, m&A- 1 ( P d i )(i),A+ 2 (pd 4 ) / 
(14) 

A drawback of this rule is that it occasionally causes motion of polymers 
that has less physical relevance. For example, when a polymer's constituent 
monomers are arranged on a straight line and all of the monomers have positive 
TKE only on the line, the polymer cannot move even if all of the monomers' 
TKE directions are identical. However, if some of the monomers have posi- 
tive TKEs in other directions, this polymer can move on average to its most 
proper direction with respect to TKE, while becoming deformed and keeping 
its integrity. 
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j=pd i {i) 



\ 



Figure 8: Translational movement and molecular bonds. If a molecule at i is 
moving to j = pd^i), (a) its bonds in directions A ±1 (pci i ) are maintained with 
their directions readjusted to A* 2 ^^), but (b) bonds in directions A ±2 (pd i ) 
and — pd i would not be maintained. 



3.4.3 Self-organization and reversibility 

In irreversible models such as the LMA, the tendency of self-organization from 
disordered high-entropy states to ordered low-entropy structures is embedded in 
their information-losing transition rules themselves. To realize "apparent" self- 
organization of ordered structures by a reversible rule without information loss, 
the rule needs additional degrees of freedom that work as a heat bath into which 
the entropy generated in the organization process should be disposed of. The 
deterministic Ising model[8] and the reversible generalization^] of the diffusion 
limited aggregation (DLA) modcl[26] are examples of this approach. We also 
adopt this approach, using the heat particle layer as the heat bath. 

When both of the neighboring sites i and j have molecules, in advance of 
the collision (12) by V'mtc, we apply site-respective TKE-heat interaction ?/>th 
defined as follows. For the molecule at site x (i.e., either i or j), when it has 
positive TKE in one direction I out of {±pd{\ (i.e., along the line connecting the 
two molecules), and when the heat particles satisfy h x j < H max and hx^i = 
at site x, then the positive TKE is transformed into a heat particle in the same 
direction I (heat release): 

mti ^ V, mtj ^ V then for x 6 {i,j} 



On the other hand, when the molecule at x does not have a TKE in either 
of {±prf i } and heat particles exist in only one of the two directions, one heat 




(15) 
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particle is transformed into the TKE in the same direction (heat absorption): 







fnti 7^ V, mtj 7^ V then for x £ {i,j} 

tke x ^± p d i \ — ( tke x> ± p d. + 1 \ if {tke x ^d il tke x ^ p d i ) = (0,0) 
h x ,± P d i J \ h x ,± p d. — 1 y and h x ^± p d. > 0, h x ^ Tp d. = 0, 

V'th (wi) = ^ otherwise. 

(16) 

Introduction of this TKE-heat interaction makes the interaction between neigh- 
boring molecules non-elastic. 

For molecular translation, potentials work as yet another control signal for 
the permutation in ipmtc of Eq. (11), and the translation is executed only if the 
changes in potentials entailed by the movement of molecule from i to j, given 
that their neighbors are fixed, can be compensated for by emission/absorption 
of heat particles at i and j: 

fnt% 7^ V, mtj = V => 

Proceed to (13) if < h u - <j(-)y*' JW < H max and 

< hjj - 5^V Sm < ff max for I £ L, 
4>mtc (w* ,Wj) = (uii , ujj ) otherwise. 

(17) 

Here, S^V*' represents the potential change in direction I caused by remov- 
ing the molecule from its current site i, and S^V''' 1 3 represents a potential 
change in direction I caused by placing the molecule (while maintaining its ori- 
entation) at vacuum site j. Further, when the translation is actually induced 
by the permutation (11), it is followed by the potential change compensation 

^pec- 

^pcc {{h u )u (h jt i)i) = (Qh,i - 6^V iAi) ) h {h jtl - ^V'^);) . (18) 

Thus, a molecule moving to a more stable site S^V*' 1 J < 0) releases 
heat particles in total, and it can be dissociated again from its neighbors only 
when enough energy is supplied by the heat particle layer. 



3.4.4 Rotation of polar molecules 

Polar molecules such as water and hydrophilic monomers can take different po- 
tential values depending on their orientations. According to the LMA rule, the 
polar molecules are rotated irreversibly into their more stable (lower-potential) 
orientations. To maintain reversibility and at the same time to enable relaxation 
into a more stable direction configuration, we utilize RKE and also the heat par- 
ticle layer. Similar to the paired site interaction in Fig. 7, rotational update is 
also performed in the interleaved scheme shown in Fig. 9. Here, we again use 
the site groups of Eq. (9) (although for the rotational permutation given below, 
interference can be prevented if the sites in each group are separated by more 
than a unit distance). 
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start 




for g in site groups G do 




for every i in a group g do 




if rati e {W, 1} then 


/ / polar molecule 


p{moi, rkei, {h,i)ieL 


) ; // rotational update 


end 





Figure 9: Interleaved update scheme for molecular orientation and RKE with 
heat interaction. 



The rotational update p is defined as follows: the polar molecule at i rotates 
according to the sign of RKE and the orientation becomes A rkei (moi) if the 
change in potentials caused by the rotation can be compensated for by emis- 
sion/absorption of heat particles at the site. Otherwise, the molecule does not 
rotate and RKE is inverted: 




5(moi, A^ 6 *)^'^); 



A rkei (m 0i 
rkei 

if < hij 
mo 

— rkei ] otherwise 
(hi,i)i 



S(moi, A r ^)^' i(i) < ffmax for V Z e L, 



(19) 



Here 



SfaA^V^W (20) 

represents the potential change in direction I that occurs when the orientation 
of the molecule at i is changed from k to A"(fc), with the neighboring molecules 
fixed. 

In this rule, the RKEs work as a kind of "second-order" signal[10], which 
preserves the history of the molecules' rotational states. Therefore, the polar 
molecules cannot stop rotating i.e., rke ^ (or given rke ^ as the initial 
condition, they cannot change their orientations forever). For an alternative 
rotation rule that allows RKEs to take values on Z, including the stationary 
state rke = 0, see Appendix A. 



3.4.5 Transportation of heat particles 

For the heat particle layer to function as a heat bath, the released heat particles 
should be effectively diffused into open areas. In our RLMA model, the diffusion 
of heat particles is conducted by a rule similar to that for the Frisch-Hasslacher- 
Pomeau lattice gas automata (FHP-LGA)[13], that is, for every unit-time up- 
date, a synchronous shift (translation) Oh is performed in each direction 

<7h : hij i ► for / e L (21) 
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followed by a local collision c/>h at each site: 



0h : 



(Oi,A'(fc))) = (m,0,0,m,0,0) 

(( h i,A'(k))) = (m,0,m,0,m,0) 
((/*i,A ! (fe))) 



i — > 



i — > 



i — > 



J (0, m, 0, 0, m, 0) if parity(pc? i ) = +1, 
1 (0, 0, to, 0, 0, to) if parity(pd^) = — 1, 

(0, TO, 0, TO, 0, TO) 

(( h i,A'(k))) otherwise 



(22) 



with < to < ifmax and fc G {1, 2, 3}. Note that the collision 0h is deterministic 
and utilizes the parity of preferential direction at each site. 

3.5 Composition of unit-time update 

After integrating the sub-steps defined above, a unit-time update of the RLMA 
can be performed through the following time sub-steps (variables in parentheses 
are those affected by the particular sub-step): 

1. Transportation of heat particles ((hi,i),ydj) 

Heat particles are diffused by the FHP-LGA-like combination of the shift 
(21) and collision (22). 

2. Rotation of polar molecules (mti, moi, rkei, (hij)) 

Using the interleaved scheme of Fig. 9, molecular rotation is performed by 
the conditional permutation (19). 

3. Molecular translation and interaction (mti, moi, (tkeij), (mbij), (fkj), pdj) 
Using the interleaved scheme of Fig. 7, 

• molecular translation is performed by the paired site conditional per- 
mutation (11) with the conditions (13) and (17), while 

• molecular interaction is performed by the conditional permutation 
(12) with the heat release (15) and absorption (16). 

4. Update of preferential direction (pdj) 

To ensure unbiasedness for the principal directions, the preferential direc- 
tion should be updated according to time. We use the simple uniform 
rotation pd i i— > A +1 (pd i ), although synchronous shifts and deterministic, 
invertiblc pseudorandom number generators can also be combined. 

These sub-steps are independent; therefore, the order can be changed. Each sub- 
time step is reversible; therefore, the inverse update is achieved by performing 
this construction in reverse. 

3.6 Conservation laws 

From the above definitions, the transition rule of RLMA conserves mass (num- 
ber of molecules) and the total energy that is given as a sum of TKEs, RKEs, 
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potential energies, and heat particles over the sites: 

i l£L i i l£L 

These conservation laws enable precise application and validation of methods 
and theorems in statistical mechanics, both equilibrium (microcanonical, canon- 
ical, and grand-canonical) and nonequilibrium (e.g., relaxation, Fourier's law of 
heat conduction, Green-Kubo relations), as is done for simpler (and in many 
cases more abstract) CA models[8, 9, 27, 28, 29, 30, 31, 32]. 

On the other hand, conservation of momenta does not hold because of 
the TKE inversion in the translational permutation (11), as mentioned in sec- 
tion 3.4.1. (Angular momenta are not conserved either, because the rotational 
permutation (19) also contains uncompensated inversion of RKE.) From the 
macroscopic viewpoint, this non-conservation of momenta seems to work posi- 
tively to enhance the model's ergodicity, instead of the effect of chaos dynamics, 
which discrete CA models lack. 

4 Simulation 

In this section, we demonstrate by simulation that the RLMA can reproduce 
the original LMA's self-organization results [3] as special cases. We also show 
that the RLMA reproduces results which arc qualitatively consistent with the 
traditional Larson-type models. More extensive results and their statistical 
treatment will be given in a future work. 

In the following simulations, the parameters are set as follows: -f/max = 8, 
Vwh-wh = %o-wo = Vwo-ip = Vip_ip = +4, Vwh-wo = Vwh-ip = —4, 

VwH-WN = KvH-IN = VwH-O = V\VO-WN = ^WO-IN = VwO-O = ^IP-WN = 

Vip-in = Vip-o = Vb-O = — 1- The simulations in section 4.1 and 4.2 adopted 
lattice space consisting of N = 24 x 24 cells, while the simulations in section 4.3 
used lattice space of N = 100 x 100 cells. Periodic boundary condition is applied 
in all the simulations, so the systems are isolated. 

4.1 Hydrophobic monomers in a polar environment 

Fig. 10 shows snapshots of the molecular layer in a simulation of a mixture of 
water and hydrophobic monomers (25% water, 25% hydrophobic monomer, 50% 
vacuum). Starting from a homogeneously mixed initial configuration with no 
heat particles (Fig. 10(a)), clustering and phase separation gradually take place 
(Fig. 10(b), (c)), accompanied by emission of heat particles (not shown in the 
figures). 

Note also that phase separation takes place in spite of the setup that the 
induction-based forces %o-h and vwo-o between a hydrophobic monomer and 
a water molecule are set equal to the dispersion interaction force woo between 
two hydrophobic monomers, and they are much weaker than the water-water 
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t = 1e+04 



t = 1e+06 



Figure 10: Snapshots of the molecular layer in a simulation of a water- 
hydrophobic monomer system. The molecules cluster preferentially with those 
of the same type, and phase separation occurs. 



binding wwh-wOi as is the case in Ref. [3] (see also section 2.3). Since the 
system is isolated, this self-organization process is entropy-driven. 

Fig. 11 shows the time evolution of mean energies per cell — TKE (sum for 
all of the principal directions QZiei ^ e i,i)»)j KKE (absolute value {\rkei\)i), 
potential energy (sum for all of the principal directions ^{J2ieL beat 
particles (sum for all of the principal directions (X^ez, — m the simulation 
run of Fig. 10. In this relaxation process, as the molecules are organized into 











Heat particles 






TKE 










RKE 






Potential 







Oe+00 2e+05 4e+05 6e+05 8e+05 1e+06 
time 

Figure 11: Time evolution of mean values of TKE, RKE, potential energy and 
heat particles per cell in a simulation of a water-hydrophobic monomer system. 

a more stable configuration, the energy released from the molecular layer is 
transferred into the heat particle layer, with the total energy conserved (mean 
total energy per cell, e to tai = E tota i/N, is 0.68). It is observed that a large part of 
the energy transfer takes place in the first few thousand steps. This corresponds 
to quick dissolution of high-potential, unstable partial configurations. 

Fig. 12 shows the time evolution of mean numbers of neighboring molecules 
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of the same types (water and hydrophobic monomer) — calculated as 

mti(i) = X}\)i s .t. mti=x for X = W, O, respectively — in the same simulation. 

The mean neighboring molecules of the same types start from the initial random 




n i i i i n 

Oe+00 2e+05 4e+05 6e+05 8e+05 1e+06 
time 



Figure 12: Time evolution of mean numbers of neighboring molecules of the 
same types (water and hydrophobic monomer) in the simulation of a water- 
hydrophobic monomer system. 

configuration (where the value ~ 0.25 x 6 = 1.5 for both water and hydrophobic 
monomers) and increase relatively slowly, taking a few million steps to reach 
the equilibrium state. This result is consistent with observations of the physico- 
chemical molecular aggregation process, where small clusters arc quickly formed, 
but as the size grows, their mobility decreases and integration into larger clusters 
requires more time. 

4.2 Amphiphilic polymers in a polar environment 

Fig. 13 shows snapshots of the molecular layer in three simulations of am- 
phiphilic tetramers (each consisting of three hydrophilic monomers plus one 
hydrophobic head monomer, see Fig. 5) in solvent water, with different settings 
for the initial distribution of heat particles. 

In the simulation of Fig. 13(a), initially, there are no heat particles, represent- 
ing the "low-temperature" condition (mean total energy per cell e to tai = 0.95). 
Although the low-temperature condition is the same as the one adopted in the 
simulation of Fig. 10, this simulation requires a longer relaxation time, because 
the polymers' mobility is lower than that of the monomers (mainly because of 
the bond maintenance condition (13)). Starting from the initial condition where 
the tetramers are homogeneously distributed, they aggregate into micelle-like 
structures, their hydrophilic heads staying in contact with water and their hy- 
drophobic tails trying to cluster. The micelle-like structures are an elementary 
example of higher-order structures [1], with emergent properties such as integrity 
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f = f = 5e+04 t = 8e+06 




t = t = 5e+04 t = 5e+06 




t = f = 5e+04 t = 2e+05 

Figure 13: Snapshots of the molecular layer in three simulations of water- 
amphiphilic polymer systems with different initial distributions of heat particles 
("temperature"): (a) low-, (b) moderate-, and (c) high-temperature conditions. 



and even lower mobility. 

Fig. 13(b) corresponds to the "moderate-temperature" condition, where the 
initial hij is given randomly from [0,1] with (/ii,z)i = 1 for I G L (mean total 
energy per cell etotai = 3.93). Compared with the low-temperature condition, 
while the sizes of the organized micelle-like structures and water aggregates 
become smaller, their motion becomes faster. 

Fig. 13(c) corresponds to the "high-temperature" condition, where the initial 
hij is given randomly from [0,4] with = 2 for I G L (mean total energy 

per cell etotai = 13.01). In this condition, the molecular motion becomes even 
faster and no distinct self-organization is observed. 

Fig. 14 shows the time evolution of mean energies per cell, and Fig. 15 
shows the time evolution of mean numbers of neighboring molecules of the same 
types (water and hydrophobic monomer) in the three abovementioned simula- 
tions with the different temperature conditions. Note the difference in the time 
scales. These results indicate the temperature dependency of the molecular pro- 
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Figure 14: Time evolution of mean values of TKE, RKE, potential energy, and 
heat particles per cell in a simulation of water-amphiphilic polymer systems 
in the different temperature conditions: (a) low-, (b) moderate-, and (c) high- 
temperature conditions. 



cess. That is, at lower temperature, polymers aggregate into larger structures; 
however, the formation process takes a longer time. On the other hand, at 
higher temperature, large structures cannot be maintained while the motion of 
polymers gets faster. This kind of temperature dependency is derived (rather 
than being presupposed) in a precise manner only from dynamical models with 
reversibility and energy conservation. Fig. 15 (especially (a)) also shows that 
the aggregation of polymers is slower than the clustering of water. 

4.3 Phase separation dynamics in ternary mixtures 

To compare in more detail the behavior of our RLMA with experimental obser- 
vations and other models (especially the Larson- type ones), we conducted sim- 
ulation of ternary mixtures of water, hydrophobic monomers, and amphiphilic 
polymers, and analyzed the phase separation dynamics with different concen- 
tration and temperature (total energy) settings. 

Theories as well as successful models have shown that the phase separation 
or domain growth dynamics generally obeys dynamic scaling[33, 15, 34, 36, 22], 
where domain structure remains statistically invariant in time under rescaling 
by the characteristic length scale L, and L grows as a function of time following 
the asymptotic power law, L(t) ~ t x l z . The theories typically suggest z = 1/3, 
though the value can differ depending on the stages of phase separation. 

To check if the RLMA realizes the dynamic scaling behavior, we investigated 
the evolution of mean cluster radius. A molecule of type X (in this case, X can 
be water W, hydrophobic monomer O, or amphiphile A) belongs to a cluster 
of type X if any of its nearest neighbor are of the same type and are already 
counted as part of the cluster. Using the cluster distribution {nx(s)} s , where 
nx{s) is the number of type X clusters with size s, the mean cluster size \x of 
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Figure 15: Time evolution of mean numbers of neighboring molecules of the 
same types (water and hydrophobic monomer) in the simulations of water- 
amphiphilic polymer systems with the different temperature conditions: (a) 
low-, (b) moderate-, and (c) high-temperature conditions. 



type X is estimated as 

s=l / s=l 

where s ma x,x is the largest cluster of type X. In two dimensions, the mean 

1/2 

cluster radius of type X is estimated as Rx ~ Xx • 

Fig. 16 shows the evolution of averaged mean cluster radius (R) = [|(xw) + 
Xo)] 1 ^ 2 : with different concentration ratios of water hydrophobic monomers 
4>q 1 and amphiphilic polymers 4>a (the ratio of vacuum </>v = 0.4 is common), 
in a temperature setting (mean total energy etotai = 1-18). Power law behavior 
is observed for all the concentration ratios (<£w, 4>o, <f>K) = (0.288,0.288,0.024), 
(0.27,0.27,0.06), and (0.24,0.24,0.12). For each setting, the scaling exponent 
1/z is estimated as 0.213 ± 0.004, 0.199 ± 0.005, and 0.121 ± 0.004, by fitting 
the data within time region [10000,200000] into (R)(t) - t l / z . These estimated 
values, especially the former two (1/z ~ 0.2) are similar to the ones obtained 
by simpler Ising spin models for binary systems [34, 35, 36], but smaller than 
the theoretical value z = 1/3 which is also obtained in Ref. [22] by adding small 
amount of amphiphile into binary mixture, like in this simulation. The small 
values are possibly because the asymptotic late stage is not reached due to the 
small size of the system, but other possible reasons can be also suggested: (i) 
the existence of hydrodynamics, which is supposed to be absent in the dynamic 
scaling hypothesis [33] as well as the Larson- type model in Ref. [22], (ii) the 
constant-energy condition, where the temperature increases as the clustering 
proceeds and the dynamic exponent becomes smaller compared to the isother- 
mal condition[15], while the latter is usually used in the Larson-type models 
and other Monte Carlo methods, and (iii) other conservation laws, which also 
work to decrease the exponent [36]. The decrease in the growth rate of (R) 



22 



o 



o o 



L. 



O 



CM - 



1e+03 



5e+03 



2e+04 



5e+04 



2e+05 



time steps 



Figure 16: Time evolution of the averaged mean cluster radius (R), with dif- 
ferent molecular ratios, o represents (^w, 4>o, 4>a) = (0.288,0.288,0.024), A 
represents (</>w, 4>o, 4>a) = (0.27,0.27,0.06), and o represents (</>w, 4>o, 0a) = 
(0.24,0.24,0.12). The lines show power law relations with estimated scaling 
exponents 1/z for their slopes. 



accompanying the increase of concentration 4>a of amphiphile is consistent with 
the result in Ref. [22]. 

For different temperature settings we also calculated the equal-time structure 
factors Sx(k, t), which is the Fourier transform of the equal-time pair correlation 
function, defined as 



The equal-time pair correlation Cxx (r, t) of type X is calculated by drawing 
shells of radius r and r + 1 around each molecule of type X, counting the 
number of the same type molecules between the shells, and finally normalizing 
by dividing by r. 

Fig. 17 shows the equal-time water-water structure factor S\v(k, t) for a sys- 
tem with (0w, 4>Oi <t>A) = (0.288,0.288,0.024) and different temperature (total 
energy) settings. In the lower temperature settings as Fig. 17(a) and (b) (mean 




(25) 



Cxx{r,t) = (S Px {r,t) 6p x 0,t)). 



(26) 
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Figure 17: The equal-time water-water structure factor for a system with 
(<£w> 4>o, 4>a) = (0.288,0.288,0.024) and different temperature (total energy) 
settings: (a) e to tai = 0.37, (b) e to tai = 1.18, (c) e to tal = 5.63. S w (k,t) at times 
t = 2, 10, 20, 50, 100, 200 x 10 3 are shown. 



total energy etotai = 0.37 and etotai = 1.18, respectively), the structure factor 
develops a peak at nonzero wave number that grows in time and the position 
of the peak moves to lower k as t increases. The peak at k ~ 1 indicates that 
the system is approaching to a global separation into spanning networks. The 
growth rate of the peak is higher for the "warm" condition etotai = 1.18 than the 
"cool" condition e to tai = 0.37. Above a critical total energy (that corresponds 
to the demixing temperature), as Fig. 17(c) (etotai = 5.63), the structure factor 
does not show any structure. These are again in good agreement with the result 
in Ref. [22], as well as experiment [37]. 



5 Conclusion 

In this paper, we described the construction of RLMA, which simulates physico- 
chemical interaction of molecules and their self-organization process. The defi- 
nition of the model has shown how to eliminate the irreversibility in the original 
LMA using several techniques to construct reversible CA. Simulation results of 
RLMA dynamics have demonstrated that the RLMA can deal with broader sit- 
uations, with the original LMA's self-organization results as special cases. The 
results also showed that the RLMA reproduces qualitatively consistent results 
with the traditional Larson- type models. 

Although several reversible CA models have been proposed to simulate self- 
organization processes [8, 9], to our knowledge, this is the first deterministic CA 
model that simulates self-organization of higher-order structures, while satisfy- 
ing strict reversibility. 

Reversibility and conservation laws of the model enable precise application 
and validation of the various methods in both equilibrium and noncquilibrium 
statistical mechanics. Reversibility also enables rigorous tracking of the infor- 
mation flow driven by the dynamics, with no veiled sources or sinks. Therefore, 



24 



the model will be preferable in analyzing the self-organization and dynamics 
of multiple levels of structures from a information-theoretic viewpoint (e.g., 
Ref. [2]), as well as from the physically grounded viewpoint. 

This study focused on the process of molecular assembly. However, the model 
can be extended to incorporate chemical reactions and catalytic effects, by intro- 
ducing more types of molecules and setting proper values of excitation energies 
for the reactions, with their modulation in the existence of neighboring catalytic 
molecules. We are currently working on the construction of such a reversible 
and thcrmodynamically consistent model that realizes "protocells" [6] with self- 
maintenance of compartment structures, metabolism, and self-reproduction. 
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A Alternative implementation of molecular ro- 
tation 

Here, we present an alternative rotation rule, which enables the stationary state 
rke = for the polar molecules. Using the same interleaved scheme of Fig. 9, the 
rotational permutation p of (19) is replaced by the new conditional permutation 
p a defined below. 

First, consider that a polar molecule at i is rotating, i.e., rkei ^ 0. Then, p a 
maintains the rotation and changes the molecule's orientation by A sgn ( rfcei ) 3 if 
the magnitude of rkei is more than sufficient to compensate for the change in 
potential induced by the rotation; else, it inverts the direction of rotation if the 
magnitude of rkei is not large enough. If the magnitude of rkei is just sufficient 
to compensate for the potential change, p & cither executes the rotation and 
brings the molecule to the stationary state, or inverts the direction of rotation, 

3 The sign function, sgn(rfcei) = +1 if rkei > 0, —1 if rkei < 0, indicates the direction of 
rotation. 
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depending on some conditions to avoid irreversibility: 



rkei 



rkej 



A sgn(rfc ei )( mo .) 

sgn(rkei)5(moi, tfMrket^yi 
\rkei\ > 5(moi, A sgn< - rke ^)V i , or 
\rke-\ = 5(moi, A sgn< - rke ^)V i 

A(5(A s s n ( rfcei ) (moi), A s s n ( rfce 
\rk ei \ = 6(moi, A s s n ( r&e *))y i 

A8 (A s s nt - rke ^ (moi), A s s n ( rfce 
Aparity(pd i ) ^ sgn(r&e. 



if rfcej 7^ and 



)V % > 0, or 



)F l < 



moi 
— rkei 



\rke 
\rke 



if rkei and 

;| < 5(moi, A Bgnt ^ rke ^)V i , or 
;j = 8(moi, A as *( rke *>)V i 
AS(A sgn< - rke ^ (moi) , A sgn ( rfee< ^)V i 
Aparity(pd i ) = sgn(rfcej). 



< 



(27) 

Next, consider that a polar molecule at i is in the stationary state, i.e., 
rkei = 0. Then, p a starts the rotation if changing the molecule's orientation by 
one of the directions A ±:L induces a negative potential change. If rotations in 
both of the directions A^ induce negative potential changes, p a starts the rota- 
tion according to the preferential direction. On the other hand, if the molecule's 
orientation is at a local potential minimum, the molecule maintains its station- 
ary state: 



moi 
rkei 



if rkei = and 



= < 



A^imoi) 
TS(moi, A* 1 )^ 

S(moi, A ±1 )\/ i < 

AS(moi, A Tl )V i > 0, or 
8(mo i ,A +1 )V i < 

A5(mOi, A" 1 )^ < A parity^) = 

m ° l J if rke^ = and 

5(moi, A +1 )^ > A 8(mOi, A" 1 )^ > 0. 



±1, 



Here 



SfcA^V* = ^8{k 1 A n )V i < 1 



(i) 



(28) 
(29) 



represents the total potential change at i that occurs when the orientation of 
the molecule at i is changed from k to A n (k) (see the notation (20), too). 

The main point is that, in p a , change in RKE and not heat compensates for 
the change in potential. Therefore, the RKE layer can work as another energy 
storage. Recall that in the rotation rule (19), RKE works just as "second-order" 
signals to preserve reversibility; hence, it cannot change to 0. It should also be 
noted that the parity (7) of the preferential direction is utilized to avoid non- 
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uniqueness of the p a 's pre-images, which could be derived from unstable fixed 
points (stationary states at orientations of local maximum potential). 

One drawback of this alternative rotation rule is that the value of RKE is 
unbounded in principle; thus, the model is not a CA in the strict sense. In 
practice, however, due to the energy conservation (23), limitless divergence of 
RKE cannot occur unless an infinite amount of energy is injected into a finite 
region. 

The RLMA model with the alternative rotation permutation p a shows quali- 
tatively similar behavior. Fig. 18 shows snapshots of the molecular layer, Fig. 19 
shows the time evolution of mean energies per cell, and Fig. 20 shows the time 
evolution of mean numbers of neighboring molecules of the same types in a 
simulation of a water-hydrophobic monomer system, with the same initial con- 
figuration as that in section 4.1. 




t = t= 1e+04 t = 1e+06 



Figure 18: Snapshots of the molecular layer in a simulation of a water- 
hydrophobic monomer system, using the alternative rotation rule p a . Clustering 
and phase separation occur in a similar manner to those shown in Fig. 10. 
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Figure 19: Time evolution of mean values of energies per cell in a simulation 
of a water-hydrophobic monomer system, using p a . The mean absolute value 
of RKE (\rkei\)i fluctuates; this is in contrast to Fig. 11 where the value is 
constant. 
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Figure 20: Time evolution of mean numbers of neighboring molecules of the 
same types in the simulation of a water-hydrophobic monomer system, using 
p a . The evolution is similar to that shown in Fig. 12. 
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